Ab initio Study of Valley Line on a Total-Energy Surface 
for Zone-Center Distortions of Ferroelectric Perovskite Oxides BaTi0 3 and PbTi0 3 
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An ab initio structure optimization technique is newly developed to determine the valley line on 
a total-energy surface for zone-center distortions of ferroelectric perovskite oxides and is applied to 
barium titanate (BaTiOa) and lead titanate (PbTiOa). The proposed technique is an improvement 
over King-Smith and Vanderbilt's scheme [Phys. Rev. B 49, 5828 (1994)] of evaluating total energy 
as a function of the amplitude of atomic displacements. The results of numerical calculations 
show that total energy can be expressed as a fourth-order function of the amplitude of atomic 
displacements in BaTi03 but not in PbTiOa. 
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PACS numbers: 77.84.-s, 77.80.Bh, 63.70.+h, 02.60.Pn 
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The j4-BC>3 perovskite oxides, such as barium titanate 
(BaTiOa) and lead titanate (PbTiOa), are extremely im- 
portant ferroelectric materials. Because of their piezo- 
electricity below Curie temperatures (408 and 763 K, 
respectively 0]) and large dielectric constants, they are 
widely used in various applications and they have been 
extensively studied experimentally and theoretically. At 
room temperature, BaTiOa and PbTiOa have tetrago- 
nal (C\ v ) structures and spontaneous polarization. The 
tetragonal structure can be understood to be a result 
of displacive transitions from their high-temperature cu- 
bic {0\) structures, as depicted in Fig. ^ The spon- 
taneous polarization is a result of the displacements of 
cations A and B and three 2_ 's to the opposite di- 
rections. While these two classic examples of ferroelec- 
tric perovskite oxides, BaTiOa and PbTiOa, have similar 
properties, such as crystal structure and unit-cell volume 
(64.4 and 63.4 A 3 , respectively), they have significantly 
different structural parameters and properties, such as 
the ratio of lattice constants c/a (1.01 and 1.06, respec- 
tively), atomic displacements (in BaTiOa the atomic dis- 
placement of B =Ti is larger than that of A =Ba, but 
in PbTi0 3 that of A=Ph is larger than that of B=Ti), 
dielectric constants, and piezoelectric constants. Cohen 
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FIG. 1: Cubic crystal structure of perovskite oxide 
ABO3. The cubic-tetragonal displacive transition accompa- 
nies atomic displacements vZ, (t=A, B, Oi, O2, O3) indicated 
by arrows. 



used the all-electron, full-potential, linearized augmented 
plane-wave (FLAPW) method to study the atomic and 
electronic structures of BaTiOa and PbTiOa .2]. He 
pointed out that the hybridization between Ti 3d and 
O 2p is strong and essential to weaken the short-range 
repulsions and allow the ferroelectric transition in both 
BaTiOa and PbTiOa. He also pointed out that the hy- 
bridization between Pb 6s and O 2p is strong so that it 
modifies the ground state and nature of the transition 
of PbTiOa by (1) increasing the ferroelastic strain that 
couples with the ferroelectric distortions or (2) hybridiz- 
ing with the valence states, leading indirectly to changes 
in the Ti-0 interactions, whereas the interaction between 
Ba and O is almost ionic in tetragonal BaTiOa. Recently, 
using synchrotron orbital radiation high-energy X-rays, 
Kuroiwa et al. confirmed the covalent bonding between 
Pb-0 and the ionic bonding between Ba-0 [3j. 

To investigate ferroelectric perovskite oxides, it is 
useful to estimate the total-energy surface for zone- 
center distortions accurately under absolute zero tem- 
perature by ab initio calculations. King-Smith and Van- 
derbilt studied the total-energy surface using ultrasoft- 
pseudopotentials and a plane- wave basis set |4|]. From 
the full symmetric cubic perovskite structure, they gave 
the displacements v T a of atoms r (=A, B, Oi, O2, O3) 
in the Cartesian directions of a (= x, y, z) along the Tis 
soft-mode normalized eigenvectors £ Q as 
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with the scalar soft-mode amplitude u a . Under the con- 
dition that the strain components r\i (i = 1, . . . , 6; Voigt 
notation) minimize the total energy for each u a , they ex- 
pressed the total-energy surface as 
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FIG. 2: The total-energy surface for zone-center distortions 
of ferroelectric perovskite oxides are schematically illustrated 
with contour lines for two strains (a) 773 = and (b) 773 > 
in the two-dimensional subspace (v^,vf) of the atomic- 
displacement space (vf, vf , wf 1 , , v° 3 ). Thick solid lines 
are the valley lines for fixed 773 's. Dashed lines show the di- 
rection of the Pis soft-mode eigenvector at zero strain. Note 
that the direction is tangential to the valley line at vl = for 
773 = 0, but this is not the case for 773 / 0. 

where u 2 = u 2 + u 2 + u 2 z , E° is the total energy for the 
cubic structure, n is one-half of the eigenvalue of the soft 
mode, and a' and 7' are the constants determined from 
coupling constants between atomic displacements and 
strains. Although their expression veritably describes es- 
sential properties of atomic displacements coupling with 
strains, it may misestimate (a) the energy gain of the 
lattice distortion from the cubic structure to the tetrag- 
onal equilibrium structure due to a deviation of atomic 
displacements from the soft-mode £ Q direction as illus- 
trated in Fig. |21 as well as (b) the total energy particu- 
larly when the amplitude of atomic displacements u z is 
larger than the equilibrium value. Also, (c) a' and 7' in 
cq. strongly depend on coupling constants between 
atomic displacements and strains, but the coupling con- 
stants, B\ xx = d 3 E/driidu 2 for instance, are not always 
constant through out all the way from the cubic struc- 
ture to the distorted equilibrium structure, and thus they 
are difficult to determine. These issues, (a)-(c), arise 
because King-Smith and Vanderbilt's expression is re- 
stricted within the atomic displacements corresponding 
to the Ti5 soft-mode eigenvector at zero strain and also 
within the fourth-order function of u a , as in eq. (J2J. In 
this work, therefore, we redefine the amplitude of atomic 
displacements as 

«- = vV) 2 + (^) 2 + + (-° 2 ) 2 + (^ 3 ) 2 . 

(3) 

and evaluate the total energy as a function of u a under 
the condition that v T a and r\i minimize the total energy 
for each u a using an ab initio norm-conserving pseudopo- 



tential method and geometric optimization. Note that 
the translational displacement should be omitted both 
in eq. Q and in eq. @, i.e., the condition that 

«£ + i£ + i£ 1 +i£' + '£ 3 = o (4) 

should be satisfied. Although, below room tempera- 
ture, BaTiOa exhibits tetragonal-orthorhombic (C^,)- 
rhombohcdral (Cf„) structure transitions and the total 
energy can be defined as a function of u x ,u y , and u z , we 
restrict our argument, in this paper, only to the cubic- 
tetragonal distortion; E to t = E to t(u x = 0,u y = 0,u z ) 
and we optimize the strains 771 = 772 and 773 while 774, 775, 
and 776 are kept at zero. 

Minimization of total energy under the constant-w a 
constraint, i.e., on the constant-7j Q sphere, is carried out 
iteratively. Suppose that at the previous iteration in- 
dexed k we calculated a gradient g, i.e., forces exerted 
on atoms, and a symmetric Hessian matrix B, i.e., the 
second derivative of total energy B aT — d 2 E to t/dv z dv T zl 
at displacement vector v^fe. In the neighborhood of v Z} k, 
total energy has the form 

-Etot(v z ) = 2*( v * _ ^z,k)B(-v z - v z ,fc) 

+*(v z - v 2:fc )g + £tot(v z ,fc) , (5) 

where t's indicate the transpose of vectors. We first min- 
imize eq. (j2Jl on a plane that includes v z ,k and is perpen- 
dicular to n = v Z) fe/|v Z) fe|, i.e., on a plane that touches 
the constant-iiQ, sphere at v Zj fc. This constraint can be 
expressed as 

*(v 2 - v Zifc )n = . (6) 

The method of Lagrange multipliers specifies that the 
gradients of eq. JHJl and eq. JHJ must be proportional. 
Thus, a trial displacement ~v' z k+1 for the next iteration 
k + 1 can be determined by solving 

B(v z - v 2 , fe )+g = An , (7) 

where A is a Lagrange multiplier. We can determine the 
Lagrange multiplier A = t nB~ 1 g/ t nB~ 1 n, and now we 
obtain 

< fe+1 =v, )fc -^(g-^|n), (8) 

but \v' z k+1 \ may be slightly different from u z . Hence, we 
use 

v' 

z,k+l 

v 2 ,fe+i = T-p ;U Z (9) 

as a trial displacement for the next iteration. The strains 
Vi = V2 and 773, i.e., lattice constants a and c, are op- 
timized simultaneously in the iterations. We continue 
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this iterative scheme of optimizing the atomic and lat- 
tice structure until the differences in total energies be- 
come less than 10 ~~ 7 Hartree twice successively. 

In the case of chemical reactions of molecules, a re- 
action path can be defined as the steepest-descent path 
from a transition structure down to reactants and down 
to products, i.e., the valley line 0- Although there 
are some techniques to determine the chemical reaction 
path 0, we expect our optimization technique to be ac- 
curate enough and to give a simple approximation for the 
ABO3 perovskite oxides. 

We use the ABINIT package |7j for all ab initio calcu- 
lations. Bloch wave functions of electrons are expanded 
into plane waves with a cut-off energy of 60 Hartree, usin 
Teter's extended norm-conserving pseudopotentials 
The pseudopotentials include O 2s and 2p, Ti 3s, 3p, 
3d and 4s, Ba 5s, 5p and 6s, and Pb 5d, 6s and 6p as 
valence electrons. Bloch wave functions are sampled on 
an 8x8x8 grid of fc-points in the first Brillouin zone. 
The grid is reduced to 40 irreducible fc-points under C\ v 
symmetry. The exchange-correlation energy is treated 
within the local density approximation (LDA). As the 
parametrized correlation energy, we use Teter's rational 
polynomial parametrization |9j , which reproduces the re- 
sults obtained by Ceperley and Alder . The electronic 
states are calculated by the iterative scheme to reach a 
tolerance of convergence that requires differences of forces 
to be less than 5x 10~ 7 Hartree/Bohr twice successively. 

The presently calculated ratios of lattice constants 
c/a and the amplitudes of atomic displacements u z for 
the tetragonal equilibrium structures are c/a = 1.0029 
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and u 
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0.0975 Bohr for BaTi0 3 and c/a = 1.0284 
0.6087 Bohr for PbTi0 3 . The experimen- 



tally observed values are c/a = 1.0086 and u z . cq = 
0.2638 Bohr for BaTi0 3 [ll| and c/a = 1.0649 and 
Mz,eq = 0.9076 Bohr for PbTi0 3 3]. Thus, calculations 
underestimate the experimental values. Calculated en- 
ergies gained through cubic to tetragonal distortions are 
0.386 meV in BaTi0 3 and 37.5 meV in PbTi0 3 . Al- 
though energy gain and experimentally observed transi- 
tion temperature (408 K = 35.2 meV for BaTi0 3 and 
763 K — 65.8 meV for PbTi0 3 ) cannot be compared di- 
rectly, the energy gains also seem to be underestimated. 
These underestimations are because the large sensitivity 
to volume makes the volume errors in LDA unusually 
important Q- Nevertheless, we believe that calculations 
using LDA clarify some trends of displacive transitions 
of perovskite oxides. 

Figure shows the obtained results for BaTi0 3 . We 
can see that total energy is well expressed as a fourth- 
order function of the amplitude of atomic displacements 
u z (Fig.[3fa)). Lattice constants a and c are also well fit- 
ted with the quadratic function of u z (Fig. E^b)). These 
results suggest that for the compound with the small 
spontaneous distortion the direction of atomic displace- 
ments remains almost equal to the soft-mode eigenvector 
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FIG. 3: (a) Calculated total energy in meV as a function 
of u z in Bohr for BaTi03 (+'s connected with solid lines) 
compared with a fourth-order function that passes through 
the calculated energy minimum at u z = 0.0975 Bohr and 
the maximum at u z = (dashed line). Zero of the energy 
scale is placed at the total energy of the cubic structure when 
u z = 0. (b) Lattice constants a and c in Bohr as functions 
of u z fitted by quadratic functions drawn with a dotted line 
and a dashed line, respectively, (c) Atomic displacements 
"f ,t)f i"? 1 = v® 2 ,v° a as functions of u z . (d) Normalized 
atomic displacements v T z /u z as functions of u z . The Tis soft- 
mode eigenvector £ z calculated by the frozen phonon method 
is additionally shown at u z = 0. 



£ z through out the speculating range of atomic displace- 
ments, as shown in Fig. |3{d), and thus the total energy 
can be well expressed within the fourth order expansion. 

In PbTi0 3 , on the contrary, total energy cannot be 
expressed as a fourth-order function of u z , as shown in 
Fig. EI a )- We can see in Fig. 0{b) that lattice constant c 
is not well fitted by the quadratic function of u z . These 
are due to the deviation of atomic displacements from the 
soft-mode £ Q direction as u z becomes larger (Fig. ^d)); 
from u z — to around u z = 0.6087, at which the total 
energy becomes minimum, the normalized displacement 
of A —Yh increases and that of B =Ti decreases. Then, 
beyond u z — 0.6087, that of 3 decreases while that of 
Oi slightly increases. This result is the first clarification 
of the non-fourth-order behavior of the total-energy sur- 
face, though it was predicted for the compounds with the 
largest spontaneous distortion 0], but has been difficult 
to formalize with coupling between atomic displacements 
and strains. 
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FIG. 4: (a)-(d) Same as Fig. |3] except the energy minimum 
is at u z = 0.6087 for PbTiOs. Note that scale of energy and 
displacements are larger than that in Fig.|!|]and that the order 
of atomic displacements of PbTiOs in (c) and (d) is A, B, 3 , 
and Oi, from top to bottom, which is different from that of 
BaTiOs. 



The conventional Slater mode 0| £siatcr = '(0 3 —1 
— 1 — l)/\/l2, is the dominant component of atomic dis- 
placements through out the entire speculating range in 
BaTiOa, as it is well known. We find that it is useful 
to analyze normalized atomic displacements of PbTiOs 
with new three components, the A-Oi,02 closing mode 
£4- 0l ,o 2 = *(2 -1 -1 0)/V6, the B-O3 closing 
mode £b-o 3 ='(0100 — l)/V2, and the rest mode 
£ rcst ='(-2 3 -2 -2 3)/\/30. The dominant A-Oi,0 2 
closing mode has its maximum at around u z = 0.6087, 
at which the total energy becomes minimum, as shown 
in Fig[S] This can be understood as follows: from u z = 
to u z = 0.6087 the covalent bonding between A-Oi,02 
becomes stronger, but beyond u z = 0.6087, rigid-sphere- 
like ionic repulsion between them becomes stronger. It 
is consistent with Cohen's description and Kuroiwa et 
al.'s experimental work 

In summary, we contrived an ab initio structure op- 
timization technique to determine the valley line on a 
total-energy surface for zone-center distortions of ferro- 
electric perovskite oxides and applied it to BaTiOa and 
PbTiC>3. (a) Our new optimization technique picks up 
the tetragonal equilibrium structure and gives its correct 
energy gain in principle, (b) The results of numerical 



calculations show that total energy can be expressed as 
a fourth-order function of the amplitude of atomic dis- 
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FIG. 5: Normalized atomic displacements of PbTiOs are de- 
composed into three modes: the j4-Oi,02 closing mode (solid 
line), the B-O3 closing mode (dashed line), and the rest mode 
(dotted line) as functions of u z in Bohr. 



placements u z in BaTi03 but not in PbTiC>3. (c) Our 
new optimization technique can automatically evaluate 
the total energy as a function of u z and coupling between 
atomic displacements and strains, i.e., the quadratic or 
nearly quadratic u Q -dependences of strains. This is an 
advantage of our technique compared with King-Smith 
and Vanderbilt's scheme. 

Computational resources were provided by the Center 
for Computational Materials Science, Institute for Mate- 
rials Research, Tohoku University. 



[1] 
[2] 
[3] 



[7] 



[8] 
[9] 

[10] 

[11] 
[12] 



E. C. Subbarao, Ferroelectrics 5, 267 (1973). 

R. E. Cohen, Nature 358, 136 (1992). 

Y. Kuroiwa, S. Aoyagi, A. Sawada, J. Harada, E. Nishi- 

bori, M. Takata, and M. Sakata, Phys. Rev. Lett. 87, 

217601 (2001). 

R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 49, 
5828 (1994). 

K. Fukui, Acc. Chem. Res. 14, 363 (1981). 
H. Bernhard Schlegel, in Modern Electronic Structure 
Theory, Part I, edited by D. R. Yarkony (World Scien- 
tific, Singapore, 1995), chap. 8. 

X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, 
M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, 
G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, 
P. Ghosez, J.-Y. Raty, and D. C. Allan, Comput. Mater. 
Sci. 25, 478 (2002). 

M. Teter, Phys. Rev. B 48, 5031 (1993). 

S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 

1703 (1996). 

D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 
(1980). 

G. H. Kwei, A. C. Lawson, S. J. L. Billinge, and S. W. 
Cheong, J. Phys. Chem. 97, 2368 (1993). 
J. C. Slater, Phys. Rev. 78, 748 (1950). 



